Sensitivity effects of void density and arrangement in a REBO high explosive 
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The shock response of two-dimensional model high explosive crystals with various arrangements 
of circular voids is explored. We simulate a piston impact using molecular dynamics simulations 
with a Reactive Empirical Bond Order (REBO) model potential for a sub-micron, sub-ns exothermic 
reaction in a diatomic molecular solid. In square lattices of voids (of equal size), reducing the size of 
the voids or increasing the porosity while holding the other parameter fixed causes the hotspots to 
consume the material more quickly and detonation to occur sooner and at lower piston velocities. 
The early time behavior is seen to follow a very simple ignition and growth model. The hotspots 
are seen to collectively develop a broad pressure wave (a sonic, diffuse deflagration front) that, upon 
merging with the lead shock, transforms it into a detonation. The reaction yields produced by 
triangular lattices are not significantly different. With random void arrangements, the mean time 
to detonation is 15.5% larger than with the square lattice; the standard deviation of detonation 
delays is just 5.1%. 
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I. INTRODUCTION 

Heterogeneities such as inclusions, voids, cracks, and 
other defects enhance the shock sensitivity of high ex- 
plosives by causing additional shock dissipation that cre- 
ates small regions of high temperature called hotspotsU 
Chemical reactions initiated in the hotspots emit pres- 
sure waves that merge with the lead shock and strengthen 
it, so that further hotspots are created with more vigor. 
This positive feedback is the principal mechanism of the 
shock-to-detonation transition in inhomogeneous explo- 
sivefl. While the significance of heterogeneities is well 
known, which of their characteristics are most important 
are not. In particular, since the details of the growth 
of reactions from the hotspots are not well understood, 
it is not known whether hotspots act separately or if the 
spatial arrangement of hotspots determines their efficacy 

Spherical voids are an often-studied, common defect in 
explosive^2HIl They can become hotspots upon collapsing 
under shock loading, and may also cause hotspots else- 
where by their partial reflection of the lead shock and by 
emitting further shocks upon their collapse and explo- 
sion. Inert beads produce similar effects; while they do 
not collapse violently enough to become hotspots, their 
reflected shocks are not weakened by immediately follow- 
ing rarefactions and so may more easily create hotspots 
where they collide. 

Experimental, theoretical, and numerical studies have 
sought to explain the sensitivity enhancement caused by 
voids and inert inclusions. Bourne and FielcP reported 
results from shocked two-dimensional samples of gelatin 
or an emulsion explosive that had large cylindrical voids 
introduced. They observed that voids could shield their 
downstream neighbors from the lead shock, but that they 
could also effect the collapse of their neighbors by emit- 
ting shock waves when they collapsed. Dattelbaum et al.'^ 
shocked samples of nitromethane with randomly embed- 



ded glass beads or microballoons, observing that their 
presence decreased the run distance to detonation and 
the pressure dependence of that distance. The balloons 
were found to have a greater effect than the beads, and 
small beads were in turn more effective than large beads. 

Medvedev et al.l^ conducted a theoretical analysis of 
emulsion explosives with microballoons that explained 
changes in detonation velocity with microballoon con- 
centration via an ignition and growth model with a con- 
stant mass burn rate per hotspot. Bourne and Miln^ 
experimentally and computationally considered a hexag- 
onal lattice of cylindrical voids in an emulsion explosive 
or nitromethane and observed that the reactions at the 
hotspots accelerated the shock relative to a comparison 
with water. In a previous reporlfiSlwe used molecular dy- 
namics (MD) to simulate single circular voids (and their 
periodic images) in a two-dimensional model solid explo- 
sive; the one rank of voids was able to induce a shock-to- 
detonation transition in the downstream material. In this 
study, we extend those simulations to samples with struc- 
tured and unstructured two-dimensional arrangements of 
voids. 



II. METHOD 

We simulate piston impacts on a number of samples 
with several equal-sized circular voids either randomly 
placed or in a regular square or triangular lattice. We 
parameterize the possible lattices by their symmetry, the 
total porosity p (proportion of molecules removed from 
the perfect crystal), and the radius r of each void. The 
spacing between voids in the square lattice is then 

5 = r\pK~fp. (1) 

The principal goal is to determine which of these param- 
eters have a significant effect on the sensitivity of the ex- 
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plosive, as measured by the time tu from piston impact 
to detonation transition. We also look for nonadditive 
contributions from the arrangement of voids (rather than 
their simple number density) and explore the mechanism 
of the development of detonation. 



A. Model 

The Reactive Empirical Bond Order (REBO) "AB" 
potential (originally developed irf^^tHI) describes an 
exothermic 2AB — A2 + B2 reaction in a diatomic 
molecular solid and exhibits typical detonation proper- 
ties but with a sub-micron, sub-ns reaction zone that 
is amenable to MD space and time scales. Heim et al. 
modified it to give a more molecular (and less plasmalike) 
Chapman- Jouguet state^^. We utilize the SPaSM (Scal- 
able Parallel Short-range Molecular dynamics) codd^ 
and the modified REBO potential ( "ModellV" ) also uti- 
lized in our previous studji^^l. The masses of A and B 
atoms are 12 and 14 amu; a standard leapfrog- Verlet in- 
tegrator is used with a fixed timestep of 0.509 fs in the 
NVE ensemble. 

Our two-dimensional samples are rectangles of herring- 
bone crystal with two AB molecules in each 6.19 x 4.21 A 
unit cell. The shock propagates in the +z direction; the 
samples are periodic in the transverse x direction. Each 
circular void is created by removing all dimers whose mid- 
points lie within a circle of a given radius. The atoms are 
assigned random velocities corresponding to a tempera- 
ture of 1 mK and an additional bulk velocity Vz — ~u,p di- 
rected into an infinite-mass piston formed by three frozen 
layers of unit cells at the —z end. The temperature is 
chosen to be small to avoid significant thermal expan- 
sion of the crystal but nonzero to avoid spurious effects 
from a mathematically perfect crystal; the RMS atomic 
displacement it causes is 2.1 pm. 

Each simulation is run until the shock (whether react- 
ing or not) reaches the free end of the sample; the traver- 
sal time of the shock (assuming that it does not acceler- 
ate) is tt = Z/us{up), where Us{up) is the shock velocity 
Hugoniot. It is broadly similar throughout the simula- 
tions, so the determination of whether or not a detona- 
tion occurs is meaningfully consistent. In particular, in 
each of the principal studies the sample length Z is held 
fixed so that tt is constant for each Up. Analysis of the 
results, including dynamic identification of molecules, lo- 
cation of the lead shock, and determination of detonation 
transition times, follow^l, except that when finding the 
detonation transition we keep the shock positions 15 nm 
apart for greater noise resistance. 



B. Systems 

We consider 94 square lattices: 27 with an integer num- 
ber n of voids (in each periodic image) and 67 with n 
allowed to merely approximate an integer (so the last 



void's distance to the end of the sample differs slightly 
from the first's to the piston). For brevity, we will term 
these two cases Si and S2 respectively. 

In Case Si, each combination of p e {1.0,1.778,2.25}% 
and r € {3,4,6} nm is considered, with Z = 1.28 jim 
chosen so that n = Z/5 is always an integer (n £ 
[12,36]), and each choice is simulated at each piston 
velocity Up G {1.96, 2.95, 3.93} km/s; 694776-2070048 
atoms are simulated. In Case S2, ^ = 845 nm and 
Up = 2.95 km/s are fixed, and the 67 p-r pairs from 
{1.0, 1.1, . . . , 2.0}% X {15, 16, ... , 59} A that yield an n 
within 0.075 of an integer are simulated {n £ [8,42], 
260736-1341296 atoms). The fixed velocity and loosened 
constraint on n allow this study to explore the p-r space 
more effectively. 

In an ancillary study called Case T, rectangular and 
triangular lattices of 71 = 10 voids are each simulated 
27 times with a fixed Up = 1.96 km/s and every com- 
bination of p e {1.0,1.5,2.0}%, r e {3,4, 5} nm, and 
Z E {416, 521, 627} nm. Here the void spacings are 
Sz = Z/n and 6^ = Trr^/pS^ = nnr^/pZ; 214060-1206248 
atoms are simulated. 

In the random case, a fixed sample size of 201 nm x 
1015 nm is used with three different (p, r, Up) triples 
taken from the Case Si lattices (but with more voids 
because of the increased sample area): Case Ri with 
(1%, 6 nm, 1.96 km/s) and n = 18, Case R2 with 
(2.25%, 6 nm, 1.96 km/s) and n = 40, and Case R3 with 
(2.25%, 4 nm, 2.95 km/s) and n = 90. For each set of 
parameters, 10 simulations are run with different ran- 
dom arrangements of voids chosen by the simple rejection 
method such that all void centers are at least 2r away 
from cither surface of the sample and at least 4r away 
from each other. This largest case involves ^3.1 million 
atoms. 



III. RESULTS 

The transition times for the 53 (of 67) Case S2 samples 
that detonated are shown in Fig. [1] The pla ne is a fit 
to the data; its equation is ^^(p, r)/ps = 9.718r/nm — 
1506p -I- 45.67. The cases that did not detonate be- 
fore the shock reached the free surface {tt « 78.2 ps) 
would occupy the rear corner of the plot (smallest p and 
largest r). The proportion of atoms that form product 
molecules was generally 85% but dropped to 65% in that 
corner. The half of those 14 simulations closer to the 
ones which detonated ended with detonation evidently 
imminent, but they are not counted as having detonated 
since a transition time could not be identified. 

At early times, we observe that the extent of reaction 
closely follows a very simple ignition and growth model. 
Suppose that a reacting hotspot is a disk that is created 
with a finite radius vq upon collapse of the void and grows 
at a constant speed v into the surrounding unreacted 
material until it overlaps its periodic images. If the disk 
contains a constant areal number density Ua of reacted 
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FIG. 1: Detonation transition times to from Case S2. The 
grid shows every p and every fourth r value and is a planar 
fit to the data, which are shown as crosses connected to it by 
lines to place them in space and indicate the discrepancies. 
The upper limit of the plotting region is the length of the 
longest, non-detonating simulations (78.2 ps). 
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FIG. 3: (Color online) Temperature over space and time in a 
sample with a square lattice of 42 voids: p — 2%, r = 16 A, 
Up = 2.95 km/s. Only the region of space occupied by the 
sample at its final compression is shown. Each spire at the 
lower left is a hotspot; the very hot region at large z, bounded 
below by a much faster shock, is the detonation. Note the left- 
moving shock generated at the transition and the pressure 
wave (apparent as a temporary strong advection) visible in 
the shocked material between z = 150 nm and 2 = 400 nm. 
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FIG. 2: Growth of reaction from the first two voids in a sample 
with p = 1% and r = 59 A. The values from the ignition and 
growth model applied to the first void and to the first two 
voids are also shown; the latter curve is almost everywhere 
indistinguishable from the simulation values. 

atoms, the number of reacted atoms as a function of time 
since initiation has the form 

N{t) = n,i?(t)2 = naTT{ro + vtf . (2) 

In Fig. [2] are plotted the counts of reacted atoms from 
the beginning of the least reactive Case S2 simulation, 
and the results of fitting one and two copies of N{t) to 
them, corresponding to the first and then the second pe- 
riodic line of voids being ignited. The two copies use the 
same growth parameters and are merely each shifted in 
time to match the data. Other Case S2 simulations have 
similar behavior, but the growth parameters depend in 
an unknown fashion on r and Up, so we have no general 
model for N{r, Up,t). 

Fig. [3] shows the temperature distribution history from 
the other extreme Case S2 simulation, with the smallest 
(eligible) r at the largest p. Each point in the figure 
is calculated with respect to the center of mass motion 



of, and averaged over, a column of computational cells 
of width Az « 0.53 nm. We note here the develop- 
ment of a broad pressure wave in the shocked material, 
visible both as a region of strong advection intersecting 
the detonation transition and as a temperature increase 
between the last few hotspots created before the tran- 
sition. It appears that the pressure waves emitted by 
the first several hotspots merge and the combined wave 
strengthens itself by encouraging the deflagration at each 
hotspot it encounters. When this wave overtakes the 
lead shock, its particle velocity is approximately equal to 
Up — 2.95 km/s, so the relative velocity in the collisions 
at the shock doubles and detonation begins immediately. 

All but 3 of the 27 Case Si simulations produced a 
detonation; those that did not used the lowest Up and 
(again) occupied the —p/+r corner of that slice of the 
parameter space. The transition times for the rest are 
shown in Fig. |4j they follow the pattern of Case S2 with 
the unsurprising addition that dto jdup < 0. The middle 
surface has the same Up as Case S2 and shows sonic non- 
planarity beyond the range of Fig. [T] With appropriate p 
and r, we see detonations even at Up < 2 km/s, which is 
much smaller than any value observed to trigger detona- 
tion with merely one periodic rank of voidg^ the feed- 
back is much strengthened by the subsequent hotspots. 

In Case T, the rectangular and triangular lattices did 
not differ significantly: they produced similar reaction 
yields and no detonations (presumably because they were 
relatively short; the Case Si simulations corresponding to 
the most reactive case treated here detonated after 87- 
90 ps). We would expect the triangular case to have a 
greater reactivity than the rectangular because the over- 



FIG. 4: Detonation transition times to from Case Si. The 
missing points for the smallest Up indicate failures to detonate. 



lap between hotspots in adjacent columns is delayed by 
their offset. The difference is never more than 5% of the 
area, however, and lasts only until the overlap is com- 
plete (about 30% of the reaction phase), so it may be 
difficult to measure. 

For Case R3, a consistent transition time of tjj = 
59.3 ps ± 5.1% was observed. The mean is 15.5% larger 
than the t]j from the corresponding Case Si simulation; 
the small standard deviation suggests that the details of 
the void arrangement are not significant. Furthermore, 
examination of one such simulation shows that a tight ar- 
rangement of 5 voids approximately 30% of the way down 
the sample produces no extra reactions. Later, a trian- 
gle of 6 voids triggers the transition, but spontaneous 
reactions elsewhere along the shock are also doing so si- 
multaneously, as shown in Fig. [5] None of the other ran- 
dom arrangements detonated. Case R2 produced yields 
of 68.0% ±3.7% (where the second percentage is relative) 
whereas its Case Si counterpart detonated. Case Ri pro- 
duced 52.1% ± 6.3% as compared to 70.1% for its coun- 
terpart. These differences are between normalized quan- 
tities, yet are partially due to the fact that the Case Si 
samples were 26.1% longer. 



IV. DISCUSSION 

It is to be expected that the detonation transition time 
ti) decreases with an increase in either the porosity p 
or piston velocity Wp, but that it increases with radius 
[dto/dr > 0) deserves further consideration. First it 
should be noted that {dtD/dr)s < 0: enlarging each of 
a set of voids without moving them (i.e., keeping the 
separation 6 fixed) does enhance the reactivity. However, 
when holding p fixed the reduction in number density 
overwhelms the effect of increasing the void size. 

The simple ignition and growth model used earlier pro- 
vides an explanation. Before any detonation begins, any 



FIG. 5: (Color online) Snapshot from a Case R3 simulation 
just as the upward-propagating shock is becoming a detona- 
tion. The whole width of the sample but only one ninth of its 
(original) length is shown. Blue atoms are unreacted, green 
are reacted, and red and purple are intermediate states. The 
red disk is a fresh hotspot; note the isolated reactions close 
to the shock everywhere. 



point in the material is reacted if and only if it is closer 
to the location of a void collapse than the R{t) associ- 
ated with that voicP^. Since all voids near a given point 
will collapse at nearly the same time, what matters is 
the expected minimum distance (dmin) to a void (after 
shock compression). We expect that (dmin) ex. 5 ex. 
so the material will react sooner, on average, with small 
voids — until, of course, the voids become so small that 
they no longer reliably produce any reactions at all. (For 
the smallest Up considered here, that failure radius is ap- 
proximately 2 nrrPSl.) 

A caveat is that if tq or w is a strong function of r, 
the larger average distances from larger voids might be 
ovewhelmed by more vigorous growth of the hotspots cre- 
ated by larger voids. While we do not have a model for 
ro(r. Up) and v{r,Up), they appear to be weak (perhaps 
sublinear) functions of r, in which case the conclusion 
of more reactivity from smaller voids holds. That v is 
a function of r at all is interesting; we suppose that the 
r-dependent strength of the reshock emitted when the 
void collapses and explodes in place may imprint on its 
surrounds a memory of the void's size. 

The model also explains how disorder in the arrange- 
ment of voids increases to- Whenever, in the random 
placement of voids, two or more are placed much closer 
than 5 to one another, their hotspots will overlap very 
quickly and the total burn front area will then be reduced; 
equivalently, (dmin) is larger for a random arrangement 
than for a lattice (especially a hexagonal one) at the same 
p and r. 

The broad pressure wave created by the lead hotspots 
seems to be the principal mechanism for the detonation 
transition in this system. Its development, the identical 
growth of the first two hotspots, the similarity of the re- 
sults from rectangular and triangular lattices, the consis- 
tency among the results from random void arrangements, 
and the apparent irrelevance of void clusters all suggest 
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that the development of a detonation is a collective ef- 
fect that depends on p, r, and the regularity of the void 
arrangement but not on the details of that arrangement. 
This collectivity affords a major simplification in predict- 
ing the behavior of collections of voids: a model might 
need only r, p, and as. 

Finally, we note that the function v{r, Up), since it will 
likely dominate ro(r, Wp) and appears to depend on both 
its arguments but not on time, may prove useful as a mea- 
sure of the strength or activity of a hotspot that might be 
incorporated into an analytical reaction rate model. For 
voids of non-uniform size, it might be sufficient to con- 
sider the variation of v in calculating a point's expected 



burn time. 
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